Anaylsis of ATAC-seq in embryos depleted for Dorsal
Kaelan (1028-02098-001-001)
September 27, 2022
1. Introduction
The goal of this analysis is to compare ATAC-seq time course data generated in wild-type, OregonR embryos to the same experiments performed in embryos depleted for Dorsal (gd7). These flies are generated from a gd7/winscy, P{hs-hid}5 stock that was heat-shocked at the larval stage (see methods). We will use these data to investigate the influence of DV patterning TFs, namely Dorsal, on chromatin accessibility, and will look at this in a pattern-specific way. We will also look at ATAC-seq and MNase-seq data comparisons between gd7 and toll10b (high, uniform Dorsal; dorsal ectoderm) mutants. Note that this analysis generates Figure 5 a-c, f-g.
2. Computational Setup
#Standard packages
library(rtracklayer) ; library(GenomicRanges); library(magrittr) ; library(Biostrings)
library(ggplot2) ; library(reshape2); library(plyranges); library(Rsamtools); library(ggpubr)
library(dplyr); library(data.table); library(patchwork); library(readr); library(testit); library(cowplot)
#KNITR Options
setwd("/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/analysis/")
options(knitr.figure_dir="figures/15_dorsal_depletion_atac/", java.parameters = "- Xmx6g")
#Lab sources
source("scripts/r/granges_common.r")
source("scripts/r/metapeak_common.r")
source("scripts/r/knitr_common.r")
source("scripts/r/caching.r")
source("scripts/r/metapeak_functions.r")
#Specific sources
library(ggseqlogo)
library(BSgenome.Dmelanogaster.UCSC.dm6)
library(rhdf5)
source("scripts/r/motif_functions.r")
#Pre-existing variables
modisco_dir <- 'bpnet/modisco/fold1/'
tasks <- c('Zld','Dl','Twi','Bcd','Cad','GAF')
threads <- 5
bsgenome<-BSgenome.Dmelanogaster.UCSC.dm63. Import relevant data for this analysis
# motifs
motifs.df <- readr::read_tsv('/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/analysis/tsv/mapped_motifs/all_instances_curated_1based.tsv.gz') %>%
as.data.frame()
# ATAC replicate coverage list: for DESeq2, we will use 3 replicates of all 4 time points, both in wt and gd7 embryos.
atac.bw.gd7.list <- list(wt_1to15_rep1 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_1to15_atac_4_cutsites.bw',
wt_1to15_rep2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_1to15_atac_6_cutsites.bw',
wt_1to15_rep3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_1to15_atac_7_cutsites.bw',
gd7_1to15_rep1 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_1to15_atac_3_cutsites.bw',
gd7_1to15_rep2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_1to15_atac_5_cutsites.bw',
gd7_1to15_rep3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_1to15_atac_6_cutsites.bw',
wt_15to2_rep1 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_15to2_atac_4_cutsites.bw',
wt_15to2_rep2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_15to2_atac_6_cutsites.bw',
wt_15to2_rep3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_15to2_atac_7_cutsites.bw',
gd7_15to2_rep1 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_15to2_atac_3_cutsites.bw',
gd7_15to2_rep2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_15to2_atac_5_cutsites.bw',
gd7_15to2_rep3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_15to2_atac_6_cutsites.bw',
wt_2to25_rep1 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_2to25_atac_4_cutsites.bw',
wt_2to25_rep2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_2to25_atac_6_cutsites.bw',
wt_2to25_rep3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_2to25_atac_7_cutsites.bw',
gd7_2to25_rep1 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_2to25_atac_3_cutsites.bw',
gd7_2to25_rep2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_2to25_atac_5_cutsites.bw',
gd7_2to25_rep3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_2to25_atac_6_cutsites.bw',
wt_25to3_rep1 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_25to3_atac_4_cutsites.bw',
wt_25to3_rep2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_25to3_atac_6_cutsites.bw',
wt_25to3_rep3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/orer_25to3_atac_7_cutsites.bw',
gd7_25to3_rep1 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_25to3_atac_3_cutsites.bw',
gd7_25to3_rep2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_25to3_atac_5_cutsites.bw',
gd7_25to3_rep3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/individual/gd7_25to3_atac_6_cutsites.bw')4. Determine differential accessibility using DESeq2
4.1. Format and calculate counts using cutsite information
# Define peaks--these are non-overlapping and evenly-sized
peaks.gr <- rtracklayer::import('/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/narrowpeak/atac_orer_all_peaks.narrowpeak') %>%
GenomicRanges::reduce(ignore.strand = T) %>%
GenomicRanges::resize(2114, 'center') %>%
plyranges::mutate(peak_id = 1:length(.))
# Calculate the counts across all peaks
counts.gd7.df<-mclapply(atac.bw.gd7.list, function(x){
regionSums(peaks.gr, x)
}, mc.cores = 8) %>% as.data.frame()
colnames(counts.gd7.df)<-names(atac.bw.gd7.list)
# Set conditions for the DESeq2 model
peak_condition_gd7.df<-data.frame(condition = c(rep('wt_1to15', 3), rep('gd7_1to15', 3), rep('wt_15to2', 3), rep('gd7_15to2', 3), rep('wt_2to25', 3), rep('gd7_2to25', 3), rep('wt_25to3', 3), rep('gd7_25to3', 3)))
peak_condition_gd7.df$condition<- factor(peak_condition_gd7.df$condition, levels = c("wt_1to15", "gd7_1to15", "wt_15to2","gd7_15to2", "wt_2to25","gd7_2to25", "wt_25to3","gd7_25to3"))
rownames(peak_condition_gd7.df)<-names(atac.bw.gd7.list)4.2. Set up differential model
library(DESeq2)
dds_gd7 <- DESeqDataSetFromMatrix(countData = counts.gd7.df,
colData = peak_condition_gd7.df,
design = ~ condition)
model_gd7 <- DESeq(dds_gd7)4.3. Gather time point results
t1_res.gd7.df <- results(model_gd7, contrast = c("condition","gd7_1to15","wt_1to15"), alpha = 0.05)
plotMA(t1_res.gd7.df, ylim=c(-3,3), ylab = 'log fold change log2(gd7_1to15/wt_1to15)') 
t1.ma <- ggplot(as.data.frame(t1_res.gd7.df), aes(x = baseMean, y = log2FoldChange))+
geom_hline(yintercept = 0, color = 'black', linetype = 'dotted')+
geom_point(aes(color = padj <0.05), size = .2) +
#geom_point(color = 'gray', size = .2) +
scale_y_continuous(limits = c(-3, 3), name = 'log2(gd7_1to15/wt_1to15)')+
scale_x_continuous(limits = c(0,10000), breaks = c(0,2000,4000,6000,8000,10000))+
xlab("Mean normalized ATAC-seq signal") +
ggtitle("MA plot at 1-1.5 hr") +
theme_cowplot()
t2_res.gd7.df <- results(model_gd7, contrast = c("condition","gd7_15to2","wt_15to2"), alpha = 0.05)
plotMA(t2_res.gd7.df, ylim=c(-3,3), ylab = 'log fold change log2(gd7_15to2/wt_15to2)')
t2.ma <- ggplot(as.data.frame(t2_res.gd7.df), aes(x = baseMean, y = log2FoldChange))+
geom_hline(yintercept = 0, color = 'black', linetype = 'dotted')+
geom_point(aes(color = padj <0.05), size = .2) +
#geom_point(color = 'gray', size = .2) +
scale_y_continuous(limits = c(-3, 3), name = 'log2(gd7_15to2/wt_15to2)')+
scale_x_continuous(limits = c(0,10000), breaks = c(0,2000,4000,6000,8000,10000))+
xlab("Mean normalized ATAC-seq signal") +
ggtitle("MA plot at 1.5-2 hr") +
theme_cowplot()
t3_res.gd7.df <- results(model_gd7, contrast = c("condition","gd7_2to25","wt_2to25"), alpha = 0.05)
plotMA(t3_res.gd7.df, ylim=c(-3,3), ylab = 'log fold change log2(gd7_2to25/wt_2to25)')
t3.ma <- ggplot(as.data.frame(t3_res.gd7.df), aes(x = baseMean, y = log2FoldChange))+
geom_hline(yintercept = 0, color = 'black', linetype = 'dotted')+
geom_point(aes(color = padj <0.05), size = .2) +
#geom_point(color = 'gray', size = .2) +
scale_y_continuous(limits = c(-3, 3), name = 'log2(gd7_2to25/wt_2to25)')+
scale_x_continuous(limits = c(0,10000), breaks = c(0,2000,4000,6000,8000,10000))+
xlab("Mean normalized ATAC-seq signal") +
ggtitle("MA plot at 2-2.5 hr") +
theme_cowplot()
t4_res.gd7.df <- results(model_gd7, contrast = c("condition","gd7_25to3","wt_25to3"), alpha = 0.05)
plotMA(t4_res.gd7.df, ylim=c(-3,3), ylab = 'log fold change log2(gd7_25to3/wt_25to3)')
t4.ma <- ggplot(as.data.frame(t4_res.gd7.df), aes(x = baseMean, y = log2FoldChange))+
geom_hline(yintercept = 0, color = 'black', linetype = 'dotted')+
geom_point(aes(color = padj <0.05), size = .2) +
#geom_point(color = 'gray', size = .2) +
scale_y_continuous(limits = c(-3, 3), name = 'log2(gd7_25to3/wt_25to3)')+
scale_x_continuous(limits = c(0,10000), breaks = c(0,2000,4000,6000,8000,10000))+
xlab("Mean normalized ATAC-seq signal") +
ggtitle("MA plot at 2.5-3 hr") +
theme_cowplot()
maplots <- t1.ma + t2.ma + t3.ma + t4.ma
ggsave("gd7_effect_ma_plots.pdf", plot = maplots, path = "figures/15_dorsal_depletion_atac/", width = 25, height = 20, units = "cm")
ggsave("gd7_effect_ma_plots.png", plot = maplots, path = "figures/15_dorsal_depletion_atac/", width = 25, height = 20, units = "cm")4.4. Plot results at enhancers
# Import enhancers list that contains some DV and AP enhancers that we manually parsed down from larger enhancer lists (Cusanovich, Koenecke, Zeitlinger) for ease of viewing.
dv_enhancers_subset.gr <- read.csv(file = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/analysis/bed/enhancers/dv_enhancers_subset.csv', header = TRUE) %>% makeGRangesFromDataFrame(keep.extra.columns = TRUE)
# Join ATAC peaks with enhancers
peaks.gr$enhancer_name <- NA
enhancer.ov <- findOverlaps(peaks.gr, dv_enhancers_subset.gr, ignore.strand = T)
peaks.gr$enhancer_name[enhancer.ov@from]<-dv_enhancers_subset.gr$name[enhancer.ov@to]
# Add patterning information
peaks.gr$pattern <- NA
peaks.gr$pattern[enhancer.ov@from]<-dv_enhancers_subset.gr$pattern[enhancer.ov@to]
#Join with DE results (I used the 4th timepoint because this is the peak influence of Dorsal binding, according to Li & Eisen, 2018)
peaks_w_de.df <- cbind(peaks.gr %>% as.data.frame, t4_res.gd7.df) %>%
dplyr::mutate(has_enhancer = ifelse(is.na(enhancer_name), 'no','yes'))
#Plot with DE results
names_to_show <- c("sog_shadow", "rho", "vn", "brk_primary", "vnd_3", "twi_Ozdemir", "mir-1_Biemar", "Mef2_Nguyen", "sna-S_Perry", "tld_Kirov", "zen_dist_Doyle", "dpp_Huang", "shn", "doc2", "Asph_Ozdemir")
gd7_effect_at_dv_enhancers <- ggplot(peaks_w_de.df, aes(x = baseMean, y = log2FoldChange))+
geom_hline(yintercept = 0, color = 'black', linetype = 'dotted')+
geom_point(aes(color = padj <0.05), size = .2) +
#geom_point(color = 'gray', size = .2) +
geom_text(data = peaks_w_de.df[peaks_w_de.df$enhancer_name %in% names_to_show,],aes(label = enhancer_name, color = pattern), size = 4)+
scale_y_continuous(limits = c(-1.5, 1.5), name = 'log2(gd7_25to3/wt_25to3)')+
scale_x_continuous(limits = c(0,10000), breaks = c(0,2000,4000,6000,8000,10000))+
xlab("Average of normalized ATAC counts") +
ggtitle("Gd7 effect at 2.5-3 hr at validated DV patterning enhancers") +
theme_cowplot()
gd7_effect_at_dv_enhancers
ggsave("gd7_effect_at_dv_enhancers_25to3_w_signif.pdf", plot = gd7_effect_at_dv_enhancers, path = "figures/15_dorsal_depletion_atac/", width = 25, height = 20, units = "cm")
ggsave("gd7_effect_at_dv_enhancers_25to3_w_signif.png", plot = gd7_effect_at_dv_enhancers, path = "figures/15_dorsal_depletion_atac/", width = 25, height = 20, units = "cm")
# Print padj values
peaks_w_de.df[!is.na(peaks_w_de.df$enhancer_name),][peaks_w_de.df[!is.na(peaks_w_de.df$enhancer_name),]$enhancer_name %in% names_to_show, c(7,8,13,14)]## enhancer_name pattern pvalue padj
## 359 dpp_Huang dorsal 1.914591e-01 5.030579e-01
## 2421 sna-S_Perry ventral 1.457059e-02 9.717441e-02
## 3218 mir-1_Biemar ventral 5.159847e-10 2.841718e-08
## 5340 Mef2_Nguyen ventral 8.185295e-03 6.346864e-02
## 5535 shn dorsal 7.561317e-02 2.947899e-01
## 6265 Asph_Ozdemir ventral 6.323338e-01 8.605463e-01
## 7290 twi_Ozdemir ventral 3.741497e-04 5.336117e-03
## 7820 rho neuroectoderm 2.977623e-02 1.618169e-01
## 8459 vn neuroectoderm 2.946955e-07 1.026187e-05
## 8954 doc2 dorsal 1.293688e-01 4.039924e-01
## 13162 zen_dist_Doyle dorsal 3.723554e-01 7.004633e-01
## 15882 tld_Kirov dorsal 5.456007e-01 8.143736e-01
## 17269 vnd_3 neuroectoderm 6.379059e-02 2.634194e-01
## 18307 brk_primary neuroectoderm 5.212455e-04 7.101970e-03
## 19654 sog_shadow neuroectoderm 4.144157e-03 3.706364e-024.5. Map DESeq2 results back to peaks
# overlap with peaks
peaks_w_de_t1.df <- cbind(peaks.gr %>% as.data.frame, t1_res.gd7.df) %>%
dplyr::mutate(has_enhancer = ifelse(is.na(enhancer_name), 'no','yes'))
peaks_w_de_t2.df <- cbind(peaks.gr %>% as.data.frame, t2_res.gd7.df) %>%
dplyr::mutate(has_enhancer = ifelse(is.na(enhancer_name), 'no','yes'))
peaks_w_de_t3.df <- cbind(peaks.gr %>% as.data.frame, t3_res.gd7.df) %>%
dplyr::mutate(has_enhancer = ifelse(is.na(enhancer_name), 'no','yes'))
peaks_w_de_t4.df <- cbind(peaks.gr %>% as.data.frame, t4_res.gd7.df) %>%
dplyr::mutate(has_enhancer = ifelse(is.na(enhancer_name), 'no','yes'))
# make GRanges for export
peaks_w_de_t1.gr <- peaks_w_de_t1.df %>%
dplyr::mutate(timepoint = '1to15') %>%
makeGRangesFromDataFrame(keep.extra.columns = TRUE)
peaks_w_de_t2.gr <- peaks_w_de_t2.df %>%
dplyr::mutate(timepoint = '15to2') %>%
makeGRangesFromDataFrame(keep.extra.columns = TRUE)
peaks_w_de_t3.gr <- peaks_w_de_t3.df %>%
dplyr::mutate(timepoint = '2to25') %>%
makeGRangesFromDataFrame(keep.extra.columns = TRUE)
peaks_w_de_t4.gr <- peaks_w_de_t4.df %>%
dplyr::mutate(timepoint = '25to3') %>%
makeGRangesFromDataFrame(keep.extra.columns = TRUE)4.6. Save DESeq2 results
deseq_at_peaks.list <- list('gd7_1to15/wt_1to15' = peaks_w_de_t1.gr,
'gd7_15to2/wt_15to2' = peaks_w_de_t2.gr,
'gd7_2to25/wt_2to25' = peaks_w_de_t3.gr,
'gd7_25to3/wt_25to3' = peaks_w_de_t4.gr)
saveRDS(deseq_at_peaks.list, '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/deseq2/gd7/deseq_timepoints_at_orer_1to3_atac_peaks.RData')4.7. Map log2FC values onto motifs
# Summarize from all time points into one total de onject
peaks_w_de_t1.df$t1_log2FC <- peaks_w_de_t1.df$log2FoldChange
peaks_w_de_t2.df$t2_log2FC <- peaks_w_de_t2.df$log2FoldChange
peaks_w_de_t3.df$t3_log2FC <- peaks_w_de_t3.df$log2FoldChange
peaks_w_de_t4.df$t4_log2FC <- peaks_w_de_t4.df$log2FoldChange
peaks_w_all_time_de.df <- peaks_w_de_t4.df[,1:8]
peaks_w_all_time_de.df$t1_log2FC <- peaks_w_de_t1.df$log2FoldChange
peaks_w_all_time_de.df$t2_log2FC <- peaks_w_de_t2.df$log2FoldChange
peaks_w_all_time_de.df$t3_log2FC <- peaks_w_de_t3.df$log2FoldChange
peaks_w_all_time_de.df$t4_log2FC <- peaks_w_de_t4.df$log2FoldChange
# Map peaks with DESeq information back to motifs
motifs_vs_peaks.ov<-findOverlaps(motifs.df %>% makeGRangesFromDataFrame(), peaks.gr, ignore.strand = T)
motifs.df$peak_id<-NA
motifs.df$peak_id[motifs_vs_peaks.ov@from]<-peaks.gr$peak_id[motifs_vs_peaks.ov@to]
motifs_across_peaks.df<-motifs.df %>%
dplyr::filter(!is.na(peak_id)) %>%
dplyr::left_join(., peaks_w_all_time_de.df %>% dplyr::select(-c(seqnames, start, end, width, strand)), by = 'peak_id') %>%
as.data.frame()5. Compare wt and gd7 ATAC-seq time course data, along with the zld- ATAC, at the patterned enhancers
5.1. Import relevant data
# Import enhancers from "Genome-wide identification of Drosophila dorso-ventral enhancers by differential histone acetylation analysis" by Koenecke et al., 2016
dv_enhancers.gr <- readRDS('/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/analysis/bed/enhancers/koenecke2016_dv_enhancers_denovo_dm6.rds')
# Create coverage bigwig list
atac.total.list <- list(orer_1to1.5 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/orer_1to15_atac_combined_normalized.bw',
orer_1.5to2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/orer_15to2_atac_combined_normalized.bw',
orer_2to2.5 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/orer_2to25_atac_combined_normalized.bw',
orer_2.5to3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/orer_25to3_atac_combined_normalized.bw',
gd7_1to1.5 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/gd7_1to15_atac_combined_normalized.bw',
gd7_1.5to2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/gd7_15to2_atac_combined_normalized.bw',
gd7_2to2.5 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/gd7_2to25_atac_combined_normalized.bw',
gd7_2.5to3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/gd7_25to3_atac_combined_normalized.bw',
zldrnai_1to1.5 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/zldrnai_1to15_atac_combined_normalized.bw',
zldrnai_1.5to2 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/zldrnai_15to2_atac_combined_normalized.bw',
zldrnai_2to2.5 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/zldrnai_2to25_atac_combined_normalized.bw',
zldrnai_2.5to3 = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/zldrnai_25to3_atac_combined_normalized.bw')5.2. Calculate ATAC-seq at mesodermal and dorsal ectodermal enhancers
atac_across_denovo.df<-mclapply(names(atac.total.list), function(x){
regionSums(dv_enhancers.gr %>% resize(1000, 'center'), atac.total.list[[x]])
}, mc.cores = 4) %>% as.data.frame()
names(atac_across_denovo.df)<-names(atac.total.list)
dv_enhancers_with_all_atac.df <- cbind(dv_enhancers.gr %>% as.data.frame, atac_across_denovo.df) %>%
dplyr::select(differential_k27ac, orer_1to1.5, orer_1.5to2, orer_2to2.5, orer_2.5to3, gd7_1to1.5, gd7_1.5to2, gd7_2to2.5,
gd7_2.5to3, zldrnai_1to1.5, zldrnai_1.5to2, zldrnai_2to2.5, zldrnai_2.5to3) %>%
as.data.table %>%
melt.data.table(id.vars = c('differential_k27ac'), variable.name = 'timepoint', value.name = 'ATAC_signal')5.3. Plot boxplots for mesodermal vs. dorsal ectodermal enhancers
box_order <- c("orer_1to1.5", "zldrnai_1to1.5", "gd7_1to1.5", "orer_1.5to2", "zldrnai_1.5to2", "gd7_1.5to2", "orer_2to2.5", "zldrnai_2to2.5", "gd7_2to2.5", "orer_2.5to3", "zldrnai_2.5to3", "gd7_2.5to3")
my_comparisons <- list(c("orer_1to1.5", "gd7_1to1.5"), c("orer_1to1.5", "zldrnai_1to1.5"),
c("orer_1.5to2", "gd7_1.5to2"), c("orer_1.5to2", "zldrnai_1.5to2"),
c("orer_2to2.5", "gd7_2to2.5"), c("orer_2to2.5", "zldrnai_2to2.5"),
c("orer_2.5to3", "gd7_2.5to3"), c("orer_2.5to3", "zldrnai_2.5to3"))
dv_enhancers_with_all_atac.df$timepoint <- factor(dv_enhancers_with_all_atac.df$timepoint, levels = box_order)
# Plot each pattern separately
meso.df <- dv_enhancers_with_all_atac.df[dv_enhancers_with_all_atac.df$differential_k27ac == "Higher in Toll10b",]
dee.df <- dv_enhancers_with_all_atac.df[dv_enhancers_with_all_atac.df$differential_k27ac == "Higher in gd7",]
denovo_meso.plot <- ggplot(meso.df, aes(x=timepoint, y=log2(ATAC_signal))) +
geom_boxplot(aes(fill = timepoint)) +
stat_compare_means(comparisons = my_comparisons, method = "wilcox.test", label = "p.signif", paired = TRUE, label.y = 13)+
ylim(8,14) +
scale_fill_viridis_d()+
ggtitle("Timecourse ATAC-seq at at mesodermal enhancers") +
theme_cowplot() +
theme(legend.position="none")
denovo_dee.plot <- ggplot(dee.df, aes(x=timepoint, y=log2(ATAC_signal))) +
geom_boxplot(aes(fill = timepoint)) +
stat_compare_means(comparisons = my_comparisons, method = "wilcox.test", label = "p.signif", paired = TRUE, label.y = 13)+
ylim(8,14) +
scale_fill_viridis_d()+
ggtitle("Timecourse ATAC-seq at dorsal ectodermal enhancers") +
theme_cowplot() +
theme(legend.position="none")
denovo_boxes <- denovo_dee.plot + denovo_meso.plot +
plot_layout(widths = c(2, 2))
denovo_boxes
ggsave("atac_regionsums_at_enhancers_1000window.pdf", plot = denovo_boxes, path = "figures/15_dorsal_depletion_atac/", width = 20, height = 10)
ggsave("atac_regionsums_at_enhancers_1000window.png", plot = denovo_boxes, path = "figures/15_dorsal_depletion_atac/", width = 20, height = 10)5.4. Plot ATAC-seq at more limited list with DEE vs. NEE vs. ME information included
# Bring in enhancers
dv_w_nee.gr <- read.csv(file = '/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/analysis/bed/enhancers/dv_enhancers_w_nee.csv', header = TRUE) %>% makeGRangesFromDataFrame(keep.extra.columns = TRUE)
# Calculate
atac_across_nee.df<-mclapply(names(atac.total.list), function(x){
regionSums(dv_w_nee.gr %>% resize(1000, 'center'), atac.total.list[[x]])
}, mc.cores = 4) %>% as.data.frame()
names(atac_across_nee.df)<-names(atac.total.list)
dv_enhancers_nee_with_all_atac.df <- cbind(dv_w_nee.gr %>% as.data.frame, atac_across_nee.df) %>%
dplyr::select(pattern, orer_1to1.5, orer_1.5to2, orer_2to2.5, orer_2.5to3, gd7_1to1.5, gd7_1.5to2, gd7_2to2.5,
gd7_2.5to3, zldrnai_1to1.5, zldrnai_1.5to2, zldrnai_2to2.5, zldrnai_2.5to3) %>%
as.data.table %>%
melt.data.table(id.vars = c('pattern'), variable.name = 'timepoint', value.name = 'ATAC_signal')
# Plot
dv_enhancers_nee_with_all_atac.df$timepoint <- factor(dv_enhancers_nee_with_all_atac.df$timepoint, levels = box_order)
ventral.df <- dv_enhancers_nee_with_all_atac.df[dv_enhancers_nee_with_all_atac.df$pattern == "ventral",]
dorsal.df <- dv_enhancers_nee_with_all_atac.df[dv_enhancers_nee_with_all_atac.df$pattern == "dorsal",]
neuroectoderm.df <- dv_enhancers_nee_with_all_atac.df[dv_enhancers_nee_with_all_atac.df$pattern == "neuroectoderm",]
ventral.plot <- ggplot(ventral.df, aes(x=timepoint, y=log2(ATAC_signal))) +
geom_boxplot(aes(fill = timepoint)) +
stat_compare_means(comparisons = my_comparisons, method = "wilcox.test", label = "p.signif", paired = TRUE, label.y = 13)+
ylim(8,14) +
scale_fill_viridis_d()+
ggtitle("Timecourse ATAC-seq at mesodermal enhancers") +
theme_cowplot() +
theme(legend.position="none")
dorsal.plot <- ggplot(dorsal.df, aes(x=timepoint, y=log2(ATAC_signal))) +
geom_boxplot(aes(fill = timepoint)) +
stat_compare_means(comparisons = my_comparisons, method = "wilcox.test", label = "p.signif", paired = TRUE, label.y = 13)+
ylim(8,14) +
scale_fill_viridis_d()+
ggtitle("Timecourse ATAC-seq at dorsel ectoderm enhancers") +
theme_cowplot() +
theme(legend.position="none")
neuroectoderm.plot <- ggplot(neuroectoderm.df, aes(x=timepoint, y=log2(ATAC_signal))) +
geom_boxplot(aes(fill = timepoint)) +
stat_compare_means(comparisons = my_comparisons, method = "wilcox.test", label = "p.signif", paired = TRUE, label.y = 13)+
ylim(8,14) +
scale_fill_viridis_d()+
ggtitle("Timecourse ATAC-seq at neuroectoderm enhancers") +
theme_cowplot() +
theme(legend.position="none")
patterning_boxes <- dorsal.plot + neuroectoderm.plot + ventral.plot +
plot_layout(widths = c(2, 2, 2))
patterning_boxes
ggsave("atac_regionsums_at_patterned_enhancers_1000window.pdf", plot = patterning_boxes, path = "figures/15_dorsal_depletion_atac/", width = 20, height = 10)
ggsave("atac_regionsums_at_patterned_enhancers_1000window.png", plot = patterning_boxes, path = "figures/15_dorsal_depletion_atac/", width = 20, height = 10)6. Compare ATAC-seq and MNase-seq between dorsalized (gd7) and ventralized (toll10b) embryos
6.1. Import data necessary for this analysis
mnase.bw.list <- list(gd7 = "/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/mnase/combined/gd7_2to4_mono_mnase_combined.bw",
toll10b = "/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/mnase/combined/toll10b_2to4_mono_mnase_combined.bw")
atac.bw.list <- list(gd7 = "/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/gd7_2to4_atac_combined.bw",
toll10b = "/l/Zeitlinger/ZeitlingerLab/Manuscripts/Zelda_and_Nucleosomes/Analysis/data/bw/dm6/atac/combined/toll10b_2to4_atac_combined.bw")6.2. Prepare functions
source("scripts/r/heatmaps/heatmaps.r")
source("scripts/r/heatmaps/metapeaks.r")
source("scripts/r/heatmaps/granges_common.r")
source("scripts/r/granges_common.r")
# Functions
read_matrix <- function(gr, cov, reverse_reads=FALSE, df=FALSE, nu=25) {
if(class(cov) == "character") {
#cov <- import.bw(cov, which=gr, as="RleList") #deprecated command, rtracklayer will recognize values
cov <- rtracklayer::import(cov, which=gr, as="RleList")
}
transform_function <- if(reverse_reads) { rev } else { identity }
gr <- granges(gr) # remove metadata to speed up RleList data retrieval in Views()
o <- order(gr)
gr <- gr[o]
if(df==TRUE){
reads.list <- regionApply_list(gr,cov, as.numeric )
reads.list <-lapply(reads.list, function(reads) { approx(reads, n=nu)$y })
reads.m <- matrix(as.numeric(unlist(reads.list, use.name=FALSE)), nrow=length(reads.list), byrow=T)
if(reverse_reads == T)reads.m <- reads.m[, ncol(reads.m):1]
}else{
reads.list <- regionApply(gr, cov, as.numeric)
#rl <- as(gr, "GRangesList")
rl <- split(gr, seqnames(gr))
rl <- rl[which(lapply(rl, function(x)length(x))!=0)] #eliminate chromosomes that are not used
view <- RleViewsList(rleList=cov[names(rl)], rangesList=rl)
reads.list <- viewApply(view, function(x) { transform_function(as.numeric(x)) })
reads.m <- matrix(unlist(sapply(reads.list, as.numeric)), nrow=length(gr), byrow=TRUE)
}
reads.m[o, ] <- reads.m
reads.m
}
standard_metapeak_matrix <- function(regions.gr, sample.cov, upstream=100, downstream=100) {
regions.gr <- resize(regions.gr, width=downstream)
regions.gr <- resize(regions.gr, width=upstream + width(regions.gr), fix="end")
failed_resize <- which(width(regions.gr) != upstream + downstream)
if(length(failed_resize) > 0) {
stop("The following regions.gr elements cannot be resized due to chromosome boundaries: ", paste0(failed_resize, collapse=", "))
}
regions.p <- regions.gr[strand(regions.gr) == "+" | strand(regions.gr) == "*"]
regions.n <- regions.gr[strand(regions.gr) == "-"]
reads <- NULL
if(class(sample.cov) == "character") {
sample.cov <- rtracklayer::import(sample.cov, which=regions.gr, as="RleList")
}
if(length(regions.p) > 0) reads <- rbind(reads, read_matrix(regions.p, sample.cov))
if(length(regions.n) > 0) reads <- rbind(reads, read_matrix(regions.n, sample.cov, reverse_reads=TRUE))
stopifnot(nrow(reads) == length(regions.gr))
reads
}
total_signal <- function(cov) {
if(class(cov) == "character") {
cache.file <- paste0(cov, ".ts.rds")
if(file.exists(cache.file)) {
return(readRDS(cache.file))
} else {
cov <- check_coverage_argument(cov)
ts <- sum(as.numeric(sapply(cov, function(x) sum(as.numeric(x)))))
#saveRDS(ts, file=cache.file)
return(ts)
}
} else {
return(sum(as.numeric(sapply(cov, function(x) sum(as.numeric(x))))))
}
}
apply_order <- function(m, o) {
m[o, ]
}
prenormalize_matrix <- function(m, bw) {
m <- m / total_signal(bw) * 150 * 15 * 1e6
}
atac_tissue_heatmap <- function(gr1, bigwigs) {
gr1 %<>% resize(1, fix="center")
gd.m1 <- standard_metapeak_matrix(gr1, bigwigs$gd7, upstream=1000, downstream=1000) %>% prenormalize_matrix(bigwigs$gd7)
toll.m1 <- standard_metapeak_matrix(gr1, bigwigs$toll10b, upstream=1000, downstream=1000) %>% prenormalize_matrix(bigwigs$toll10b)
diff.m1 <- toll.m1 - gd.m1
o1 <- order(rowSums(diff.m1[, 750:1250]))
gd.m <- apply_order(gd.m1, o1)
toll.m <- apply_order(toll.m1, o1)
diff.m <- apply_order(diff.m1, o1)
list(gd=gd.m, toll=toll.m, diff=diff.m1)
}
mnase_tissue_heatmap <- function(gr1, bigwigs) {
gr1 %<>% resize(1, fix="center")
gd.m1 <- standard_metapeak_matrix(gr1, bigwigs$gd7, upstream=1000, downstream=1000) %>% prenormalize_matrix(bigwigs$gd7)
toll.m1 <- standard_metapeak_matrix(gr1, bigwigs$toll10b, upstream=1000, downstream=1000) %>% prenormalize_matrix(bigwigs$toll10b)
diff.m1 <- gd.m1 - toll.m1
o1 <- order(rowSums(diff.m1[, 750:1250]))
gd.m <- apply_order(gd.m1, o1)
toll.m <- apply_order(toll.m1, o1)
diff.m <- apply_order(diff.m1, o1)
list(gd=gd.m, toll=toll.m, diff=diff.m1)
}
# normalization in the function: you are taking the 98th percentile value (which is 47.43776; this comes from the draw_standard_heatmap function) and dividing all values in the matrix by that. You then take that value if it is less than 1 and take 1 if it is greater than 1 so that everything is technically less than 1
normalize_matrix <- function(m, value) {
m <- pmin(m / value, 1)
m
}
jeff_normalize_heatmap <- function(m, normalize=TRUE) {
if(normalize) {
max.value <- quantile(m, 0.98, na.rm=TRUE)
m <- normalize_matrix(m, max.value)
return(m)
}
}
# write function to get the order that we use, which in this case is based on ATAC
atac_order <- function(gr1, bigwigs) {
gr1 %<>% resize(1, fix="center")
gd.m1 <- standard_metapeak_matrix(gr1, bigwigs$gd7, upstream=1000, downstream=1000) %>% prenormalize_matrix(bigwigs$gd7)
toll.m1 <- standard_metapeak_matrix(gr1, bigwigs$toll10b, upstream=1000, downstream=1000) %>% prenormalize_matrix(bigwigs$toll10b)
diff.m1 <- gd.m1 - toll.m1
order_for_atac <- order(rowSums(diff.m1[, 750:1250]))
return(order_for_atac)
}6.3. Calculate and plot ATAC-seq in gd7 and toll10b embryos
# We want to use the mesodermal enhancers, which comes from the higher in toll10b Nej peaks
peaks_toll <- subset(dv_enhancers.gr, differential_k27ac == "Higher in Toll10b")
# Calculate
dv_atac_heatmaps <- atac_tissue_heatmap(peaks_toll, atac.bw.list)
gd_atac_normalized.mat <- jeff_normalize_heatmap(dv_atac_heatmaps$gd)
toll_atac_normalized.mat <- jeff_normalize_heatmap(dv_atac_heatmaps$toll)
# Format and plot
heatmap_colors3 <- colorRampPalette(c("#F6F6F6", "#330000"))(32)
gd_atac_normalized.mat.df <- gd_atac_normalized.mat %>% as.data.table
colnames(gd_atac_normalized.mat.df) <- c(-1000:(1000-1)) %>% as.character
gd_atac_normalized.mat.df <- gd_atac_normalized.mat.df %>%
dplyr::mutate(row_id = 1:nrow(.)) %>%
melt.data.table(id.vars = c('row_id'), variable.name = 'position', value.name = 'reads')
gd7_atac_at_mesodermal_regions_center_nej_peaks_norm.plot <- ggplot()+
ggrastr::geom_tile_rast(data=gd_atac_normalized.mat.df, aes(x=position, y=row_id, fill=reads))+
scale_fill_gradientn(colors = heatmap_colors3,
name = "norm.\nsignal") +
scale_y_continuous(name = 'H3K27ac higher in Tol10b (Mesodermal regions), n=416')+
xlab('Distance from Nej peak (bp)') +
#scale_x_continuous(expand = c(0, 0), name = 'Distance from Nej peak (bp)') +
ggtitle('Normalized ATAC reads in gd7 (poised enhancers)')+
theme_cowplot()+
theme(text=element_text(size=14), legend.position="bottom", panel.background = element_blank())
toll_atac_normalized.mat.df <- toll_atac_normalized.mat %>% as.data.table
colnames(toll_atac_normalized.mat.df) <- c(-1000:(1000-1)) %>% as.character
toll_atac_normalized.mat.df <- toll_atac_normalized.mat.df %>%
dplyr::mutate(row_id = 1:nrow(.)) %>%
melt.data.table(id.vars = c('row_id'), variable.name = 'position', value.name = 'reads')
tol10b_atac_at_mesodermal_regions_center_nej_peaks_norm.plot <- ggplot()+
ggrastr::geom_tile_rast(data=toll_atac_normalized.mat.df, aes(x=position, y=row_id, fill=reads))+
scale_fill_gradientn(colors = heatmap_colors3,
name = "norm.\nsignal") +
scale_y_continuous(name = 'H3K27ac higher in Tol10b (Mesodermal regions), n=416')+
xlab('Distance from Nej peak (bp)') +
#scale_x_continuous(expand = c(0, 0), name = 'Distance from Nej peak (bp)') +
ggtitle('Normalized ATAC reads in tol10b (active enhancers)')+
theme_cowplot()+
theme(text=element_text(size=14), legend.position="bottom", panel.background = element_blank())
atac_tissue.plot <- gd7_atac_at_mesodermal_regions_center_nej_peaks_norm.plot + tol10b_atac_at_mesodermal_regions_center_nej_peaks_norm.plot
atac_tissue.plot
ggsave("atac_at_mesodermal_regions_center_nej_peaks_norm.pdf", plot = atac_tissue.plot, path = "figures/15_dorsal_depletion_atac/", height = 6, width = 8)
ggsave("atac_at_mesodermal_regions_center_nej_peaks_norm.png", plot = atac_tissue.plot, path = "figures/15_dorsal_depletion_atac/", height = 6, width = 8)6.4. Calculate and plot MNase-seq in gd7 and toll10b embryos
dv_mnase_heatmaps <- mnase_tissue_heatmap(peaks_toll, mnase.bw.list)
gd_normalized.mat <- jeff_normalize_heatmap(dv_mnase_heatmaps$gd)
toll_normalized.mat <- jeff_normalize_heatmap(dv_mnase_heatmaps$toll)
# Formant and plot
heatmap_colors2 <- colorRampPalette(c("#F6F6F6", "#19334d"))(32)
gd_normalized.mat.df <- gd_normalized.mat %>% as.data.table
colnames(gd_normalized.mat.df) <- c(-1000:(1000-1)) %>% as.character
gd_normalized.mat.df <- gd_normalized.mat.df %>%
dplyr::mutate(row_id = 1:nrow(.)) %>%
melt.data.table(id.vars = c('row_id'), variable.name = 'position', value.name = 'reads')
gd7_mnase_at_mesodermal_regions_center_nej_peaks_norm.plot <- ggplot()+
ggrastr::geom_tile_rast(data=gd_normalized.mat.df, aes(x=position, y=row_id, fill=reads))+
scale_fill_gradientn(colors = heatmap_colors2,
name = "norm.\nsignal") +
scale_y_continuous(name = 'H3K27ac higher in Tol10b (Mesodermal regions), n=416')+
xlab('Distance from Nej peak (bp)') +
#scale_x_continuous(expand = c(0, 0), name = 'Distance from Nej peak (bp)') +
ggtitle('Normalized MNase reads in gd7 (poised enhancers)')+
theme_cowplot()+
theme(text=element_text(size=14), legend.position="bottom", panel.background = element_blank())
toll_normalized.mat.df <- toll_normalized.mat %>% as.data.table
colnames(toll_normalized.mat.df) <- c(-1000:(1000-1)) %>% as.character
toll_normalized.mat.df <- toll_normalized.mat.df %>%
dplyr::mutate(row_id = 1:nrow(.)) %>%
melt.data.table(id.vars = c('row_id'), variable.name = 'position', value.name = 'reads')
tol10b_mnase_at_mesodermal_regions_center_nej_peaks_norm.plot <- ggplot()+
ggrastr::geom_tile_rast(data=toll_normalized.mat.df, aes(x=position, y=row_id, fill=reads))+
scale_fill_gradientn(colors = heatmap_colors2,
name = "norm.\nsignal") +
scale_y_continuous(name = 'H3K27ac higher in Tol10b (Mesodermal regions), n=416')+
xlab('Distance from Nej peak (bp)') +
#scale_x_continuous(expand = c(0, 0), name = 'Distance from Nej peak (bp)') +
ggtitle('Normalized MNase reads in tol10b (active enhancers)')+
theme_cowplot()+
theme(text=element_text(size=14), legend.position="bottom", panel.background = element_blank())
mnase_tissue.plot <- gd7_mnase_at_mesodermal_regions_center_nej_peaks_norm.plot + tol10b_mnase_at_mesodermal_regions_center_nej_peaks_norm.plot
mnase_tissue.plot
ggsave("mnase_at_mesodermal_regions_center_nej_peaks_norm.pdf", plot = mnase_tissue.plot, path = "figures/15_dorsal_depletion_atac/", height = 6, width = 8)
ggsave("mnase_at_mesodermal_regions_center_nej_peaks_norm.png", plot = mnase_tissue.plot, path = "figures/15_dorsal_depletion_atac/", height = 6, width = 8)7. Conclusions
Here, we have characterized the differences in Ore-R and gd7 ATAC-seq and find that the changes depend on the pattern of expression. While mesodermal enhancers, which are Dorsal-activated, decrease in accessibility upon loss of Dorsal, dorsal ectodermal enhancers, which can be Dorsal-repressed, do not lose accessibility but instead can gain it. This suggests that Dorsal’s influence on chromatin accessibility is commensurate with its role as an activator or repressor of transcription. We see these differences further highlighted by comparing mesodermal enhancers in active (toll10b) and poised (gd7) states, and find that active enhancers are more accessibilie while poised enhancers have more nucleosome occupancy.
8. Session Information
For the purposes of reproducibility, the R/Bioconductor session information is printed below:
sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS: /n/apps/CentOS7/install/r-4.2.0/lib64/R/lib/libRblas.so
## LAPACK: /n/apps/CentOS7/install/r-4.2.0/lib64/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] DESeq2_1.36.0 SummarizedExperiment_1.26.1
## [3] Biobase_2.56.0 MatrixGenerics_1.8.1
## [5] matrixStats_0.62.0 rhdf5_2.40.0
## [7] BSgenome.Dmelanogaster.UCSC.dm6_1.4.1 BSgenome_1.64.0
## [9] ggseqlogo_0.1 DECIPHER_2.24.0
## [11] RSQLite_2.2.16 viridis_0.6.2
## [13] viridisLite_0.4.1 digest_0.6.29
## [15] pander_0.6.5 lattice_0.20-45
## [17] cowplot_1.1.1 testit_0.13
## [19] readr_2.1.2 patchwork_1.1.2
## [21] data.table_1.14.2 dplyr_1.0.9
## [23] ggpubr_0.4.0 Rsamtools_2.12.0
## [25] plyranges_1.16.0 reshape2_1.4.4
## [27] ggplot2_3.3.6 Biostrings_2.64.1
## [29] XVector_0.36.0 magrittr_2.0.3
## [31] rtracklayer_1.56.1 GenomicRanges_1.48.0
## [33] GenomeInfoDb_1.32.3 IRanges_2.30.1
## [35] S4Vectors_0.34.0 BiocGenerics_0.42.0
##
## loaded via a namespace (and not attached):
## [1] ggbeeswarm_0.6.0 colorspace_2.0-3 ggsignif_0.6.3
## [4] rjson_0.2.21 ellipsis_0.3.2 rstudioapi_0.14
## [7] farver_2.1.1 bit64_4.0.5 AnnotationDbi_1.58.0
## [10] fansi_1.0.3 splines_4.2.0 codetools_0.2-18
## [13] cachem_1.0.6 geneplotter_1.74.0 knitr_1.40
## [16] jsonlite_1.8.0 Cairo_1.6-0 broom_1.0.1
## [19] annotate_1.74.0 png_0.1-7 httr_1.4.4
## [22] compiler_4.2.0 backports_1.4.1 assertthat_0.2.1
## [25] Matrix_1.4-1 fastmap_1.1.0 cli_3.3.0
## [28] htmltools_0.5.3 tools_4.2.0 gtable_0.3.0
## [31] glue_1.6.2 GenomeInfoDbData_1.2.8 Rcpp_1.0.9
## [34] carData_3.0-5 jquerylib_0.1.4 vctrs_0.4.1
## [37] rhdf5filters_1.8.0 xfun_0.32 stringr_1.4.1
## [40] lifecycle_1.0.1 restfulr_0.0.15 rstatix_0.7.0
## [43] XML_3.99-0.10 zlibbioc_1.42.0 scales_1.2.1
## [46] vroom_1.5.7 ragg_1.2.2 hms_1.1.2
## [49] RColorBrewer_1.1-3 yaml_2.3.5 memoise_2.0.1
## [52] gridExtra_2.3 ggrastr_1.0.1 sass_0.4.2
## [55] stringi_1.7.8 highr_0.9 genefilter_1.78.0
## [58] BiocIO_1.6.0 BiocParallel_1.30.3 systemfonts_1.0.4
## [61] rlang_1.0.4 pkgconfig_2.0.3 bitops_1.0-7
## [64] evaluate_0.16 purrr_0.3.4 Rhdf5lib_1.18.2
## [67] labeling_0.4.2 GenomicAlignments_1.32.1 bit_4.0.4
## [70] tidyselect_1.1.2 plyr_1.8.7 R6_2.5.1
## [73] generics_0.1.3 DelayedArray_0.22.0 DBI_1.1.3
## [76] pillar_1.8.1 withr_2.5.0 survival_3.4-0
## [79] KEGGREST_1.36.3 abind_1.4-5 RCurl_1.98-1.8
## [82] tibble_3.1.8 crayon_1.5.1 car_3.1-0
## [85] utf8_1.2.2 tzdb_0.3.0 rmarkdown_2.16
## [88] locfit_1.5-9.6 grid_4.2.0 blob_1.2.3
## [91] xtable_1.8-4 tidyr_1.2.0 textshaping_0.3.6
## [94] munsell_0.5.0 beeswarm_0.4.0 vipor_0.4.5
## [97] bslib_0.4.0